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ABSTRACT 



The purpose of this thesis is to develop an algorithm for segmenting images 
corrupted by a high level of noise with different characteristics. In particular the 
images considered are composed of several regions describing different objects and 
background. The algorithm described is based on a Markov Random Field (MRfi 
model of the image and uses Kalman Filtering (KF) techniques and Dynamic 
Programming (DP) in order to smooth within the regions. The theoretical background 
for one dimensional and two dimensional data which have different characteristics and 
simulation results are presented, with examples of synthetic data and underwater 
imaees. 
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I. INTRODUCTION 



The objective of segmentation is to divide a given image into meaningful regions 
or units. In many vision problems the first task is to distinguish objects to discriminate 
between various obstacles. 

A general vision system for Artificial Intelligence applications can be decomposed 
in the following blocks: 

1. Data acquisition : The picture of interest can be taken by suitable cameras for 
these purposes and can be digitized using proper software and hardware. 

2. Segmentation : Removing noise and dividing the image into significant regions. 

3. Image understanding: Deciding which objects are present and classifying them 
according to size, shape, etc. 

The general approach at the basis of image segmentation is based on the 
identification of different properties characterizing the regions of interest. For 
example, several objects can be identified by their average intensity levels or by the 
textures on the surfaces. The task of an image segmentation algorithm is therefore to 
identify the regions from the vision signals which contain noise (due to the electronic 
equipment and other environmental factors such as turbid waters in underwater 
environments) or uninteresting details. For example, if it is desired to recognize the 
presence of a house in the picture, we need to segment the image by disregarding the 
doors, windows, cracks on the wall, etc. 

Although all the methods in the literature for segmentation are well suited to 
most of the cases of interest, they have poor performances in the presence of a high 
level of noise. In the cases of noisy data statistical techniques are more suitable. In 
particular these are based on classical techniques of estimation, such as Maximum 
Likelihood (ML) or Maximum a Posteriori (MAP). Statistical estimation techniques 
rely on the use of distinct statistical models for the original image and for the 
disturbances. 

In this thesis we present a segmentation algorithm for images affected by a high 
level of noise, based on statistical techniques. In particular we use Markov Random 
Fields (MRF) as a model of the original image, and we assume White Gaussian noise 
as a model for the disturbance. The reason for the choice of MRF is based on several 
considerations: 

1. They extend one dimensional well-known models to the multidimensional case. 
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2. They relate global statistics to local dependencies of the data. 

3. They can model compact, distinct regions by a suitable choice of parameters. 

The smoothing is obtained by a combination of MRF (for the transition between 

regions). Kalman Filtering Techniques (to filter within the regions) and Dynamic 
Programming (to determine the best sequence of edges which maximizes the likelihood 
function). 

Related to the works of segmentation, the algorithms of Derin. Elliott and Cristi 
[Refs. 1.2], Geman [Ref. 3], Besag [Ref. 4] have shown the effectiveness of these 
techniques segmenting images corrupted by a high level of noise. Parallel to these 
works, a combination of Autoregressive and MRF models has been presented by 
Therrien [Ref. 5] to segment textured images of terrain data. 

In the next chapter, the properties of modeling the original scene by Markov- 
Random Field and choice of the parameters in Markov Random Fields are presented. 
Estimation algorithms for one and two dimensional data are given in Chapter 111. 
Finally, simulation results for different characteristics of data are the subjects of 
Chapter IV. 
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II. STATEMENT OF THE PROBLEM AND MODELING USING MRF 



The problem to be addressed in this chapter is how to model the original scene 
and its estimation from the given (observed) data. As stated in the Introduction the 
model for the original scene is assumed to be a Markov Random Field (MRFi whose 
properties are discussed in the next section. As explained below, the estimation is 
based on a Bayesian approach and a Maximum a Posteriori (MAP) techniques. 

A. PROBLEM STATEMENT 

Let the original scene X = {X k T be a random Held on a Finite two dimensional 
lattice L = {(k,t) s 7 ? : 1 < k ^ Nj. 1 < t < X-,} where Z is the set of integers. 

assuming discrete values Fj, F,. F 3 which are constant in regions Rj. R,. R,.... with 

Rj c L and Rj H R. = 0 For i*j. For example this might be the case of several objects 
with different intensity levels. 

Also let X be corrupted by an additive noise and modeled by a random field 
\V={W k }, (k.t) e L. Furthermore, W will be assumed to be a White Gaussian field 
with zero mean and known variance <7 2 . W is identically independently distributed 
(i.i.d.) and independent of X. 

Therefore, the observed image can be given by a random field V = (Y k t }. as 

\.t = §k.t (\,t ’ ^ k,t^ ( -' 1 ’ 

where (k.t) e L and g k t (. , .) is an arbitrary function. In our case we assume additive 
noise, and therefore the observed signal can be expressed as 



Y„ . = X„ . + W 



' k,t 



k.t 



k,t 



( 2 .:) 



Now, the problem is to estimate X from the data Y. Since we assume that only 
noisy observations are available, the estimate can be obtained by maximization of the 
a posteriori distribution of the scene with respect to x by a Bayesian approach. 
Therefore 



P X|Y (x|y) = max (x} P X | Y (x|y) 



(2.3) 
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where x is tiie estimate of original scene and P\|y(x|y) = Probability (X = \,Y = y». 
As it can be seen front Equation 2.3. the problem is to maximize the a posters n 

s A 

probability of x given y. This form of Bayesian estimation is known as M ax. mum a 
Posteriori or MAP estimation. 

Bayes factorization and the convenience of logarithmic operation y ields 

lnPyjxOA) + Inl’y(x) = max^j [InPy^yOJx) + lnPy(x)j 1 2.-1 » 

The two likelihood terms on the left and right hand sides of Equation 2.4 are 
determined by the model of X and of the disturbance. The model for the scene X is 
assumed to be a Markov Random Field (MRF) which is formally defined by Besag 
| Ref. 3: p. 724J and by Elliot. Derm. Cristi and Geman [Refs. 0.7], 

B. MARKOV RANDOM FIELDS (MRF) 

Definition : Let L be a finite lattice L = [(k.t) : 1 < k < Nj. 1 < t < NX] and 
H k t defined as a subset of L so that (k.t) s q k t where (k.t) e L and (i.j) e . if and 
only if (k.t) 6 \\. .. Then a random field X = (X k t ] with the property 



P[X b = x. , | X. . = \. . 

1 k,t k.t 1 i.j i.j 

X. . = X. . 

i.j i.j 



P» \ = v 

1 ,A k.t x k.t 



(i.j) e Lj = 

< i *j> e 



(2.0) 



is a Markov Random Field with neighborhood q k . 

A Markov Random Field can be illustrated as in Figure 2.1, where the statistics 
of the element (k.t) indicated by a depends on its neighbors indicated by a ' O" 
only. Two simple neighborhoods are shown in Figure 2.2. These are: 



Mij 1 = UU):0 < (i-1) 2 + (j-n) 2 < 1} 



( 2 . 6 ) 



and 



= (( 1 .n) : 0 < (i-1) 2 + (j-n) 2 < 2} 



(2.7) 



Relatively simple neighborhood systems as in Figure 2.2 are adequate in modeling most 
scenes of interest. 



o o o 

0*0 

o o o 



Figure 2.1 Dependence of Markov Random Field 
P x r | Everything else). 
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Figure 2.2 Neighborhoods n. - 1 (a) and n. . 2 (b). 

Although a MRF can be considered as a multidimensional extension of a 
Markov Chain, a difficulty exists in defining MRF. In fact the main concept at the 
basis of Markov Chains is the concept of transition probabilities. These are subjects to 
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mild conditions and they rely on the fact that an ordering can be determined in one 
dimensional problems. This is not the case in multidimensional problems where an 
ordering of points in general does not exist, which presents a difficulty in the definition 
of transition probabilities. Because of this we have to approach MRFs in a different 
way from the Markov Chains. Therefore, to determine the structure for the joint 
probability of MRF models, we need some new definitions and theorems. In 
particular the joint probability is based on the concept of clique gisen below: 

Definition: Given a Markov Random field (MRF) with neighbors q a einjite is a 
set of pixels which are neighbors of each other [Ref. 6: p.l9S). The various types of 
cliques for rjj j 1 and rjj j" are shown in Figure 2.3. Based on the definition of cliques, 
the joint probability of a Markov Random Field with neighborhood lj can be expressed 
by the following theorem: 

Theorem (Hammersley and Clifford) [Ref. 6: pp. 192-199]. Suppose X is an MRF 
such that P^(x> '• 0 for all x. with neighborhood rj. Then P x (x) can be defined as 
following: 



p x (x)= Y e ' L(x) (2.S) 

where 

U(x) = £ v c (x) (2.9) 

ceg 

^ is the set of all cliques as shown in Figure 2.3 and usually L(x) is defined as the 
energy function. V (x) the potential associated with clique c. and the partition function 
Z is defined as 

z= V e -L'(x) (2.10) 

x 

which is a normalizing constant that causes P x to be a consistent probability measure. 
On the basis of the above theorem we can arbitrarily assign a joint probability of 
MRFs by the potential functions \' c (x). The only constraint on V c (.x) is that it has to 
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(a) 



□ 



| 

(b) 



Figure 2.3 Cliques for Neighborhood Systems q. (a) and q-, ■' (b). 

be Finite for all x. It means that the form of V (x) is fixed bv the structure of the 
lattice given by Kinderman [Ref. S]. 

in order to understand the above result and definitions, the joint probability of 
anv Markov Random Field with neighborhood n. 1 can be written as 

'i.j 



lnP X ( x > = Iv I (\. t ) + Iv 2( x k<l . 1 ,x k<t) + I V 3 ( x k . 1)t .x ku) -lnZ 
k.t k.t k.t 



:2.u) 



for the case in Figure 2.3 (a). 
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Definition: A MRF which has all nonzero probabilities as in liquation 2.11. 
namely, P x (x) > is called a Gibbs Field (GF). The Gibbs model can be used to 
model the spatial interactions between regions, in which a pixel and its neighbors have 
similar statistical properties. In other words, it models spatial continuity which 
characterizes the images we want to segment. Note that the .structure of P^ix) is 
relatively simple for as in Equation 2.1 1. Its complexity increases considerably for 
la i ser neighborhoods. For that reason, in General onlv the lower sets n. .’ and n ■ are 

h.j 'uj 

considered practically applicable, and even m i] i ■' cliques with only one or two 
elements are taken into account. In addition, the Gibbs distribution might also be 
used to model textures. In our problem we just want to model the fact that the most 
likely scenes have regions which are clusters of pixels, and the Markov Random Field 
model is adequate for this purpose, provided by properly assigning the potential 
functions. A particular model we are going to consider in this thesis is the MRF with 
neighbors t| l as in Figure 2.3 (a) and with potential functions defined as: 



v l< x k.i> = 0 



( 2 . 12 ) 



v : (x k.t-r\.t> = v j( x k-i.u x k.t) = 



P ifx k.t = x k,t-i 

-P otherwise 



(2.13) 



for any value of k and t and p a positive constant. The larger the parameter p, the 
more the joint distribution P x is peaked around smooth realizations. By this definition 
we can write the joint probabilities as 



lnP x (x ) = P f Hg(x kt ,x k>t . 1 ) + g(x k(t .x k . 1 t )} -InZ 
k.t 



(2.14) 



where 



S< x k.t 



k-l,i' 



1 if x. = x. , , 

k.t k-l,t 

-1 otherwise 



(2.15) 



Also, by the assumption of Gaussian noise with zero mean and standard deviation <7, 
Equation 2.14 with Equation 2.4 yields the likelihood as 
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( 116 ) 



k.t k.t 



Hence, given the data y= (y k in noise, the noise model \V(<|), <r ) and the Gibbs field 
model Equation 2.14 with parameter p, the estimate of the original scene is given by 

A 

x= (x k } such that 



In particular there is no need to undertake the difficulty of calculating the 
partition function Z since it is a normalizing constant. So. it was omitted in Equation 
2.17. Approaches to minimization of Equation 2.17 can be based on relaxation 
[Ref. 3: pp. 730-732] or deterministic [Ref. 4: pp. 266-267] or Dynamic Programming 
[Ref. 2: pp. 44-45 and Ref 8: pp. 12-13]. In the next chapter a different approach to 
minimization of Equation 2.17 based on Kalman Filtering techniques will be presented. 

C. CHOICE OF THE PARAMETER IN THE MRF 

The difficulty in using the Gibbs Distribution as a region model is to estimate the 
parameters of the model from specific realizations. Some methods have been used 
previously to estimate the model parameters. For example, Besag [Ref 6: pp. 211-212] 
suggested the coding method where the parameters are determined from subsets of data. 
This requires solution of a set of nonlinear equations. Another example of the 
parameter estimation method is given by Derin and Elliot [Ref. 2: pp. 43-44] which 
uses standard linear, least-squares estimation. This has been proved to be effective in 
modeling textures by GFs. A different approach was used by Geman [Ref 3: pp. 
727-729] which defined a simulated annihiling, where stochastic relaxation was 
combined with monotonic increasing parameters. 

For this thesis, the parameter P of the model in Equation 2.17 was set by trial 
and error until a reasonable filtering was reached. The optimum value of the 
parameter P is smaller for high values of signal to noise ratio. So, it is necessary to 
modify this parameter as the signal to noise ratio changes. 




(2.17) 



k.t 



k.t 



is minimal. 
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III. SEGMENTATION AND SMOO THING OF ONE DIMENSIONAL AND 

TWO DIMENSIONAL SIGNALS 



The content of this chapter is the segmentation algorithms lor one dimensional 
and two dimensional signals. Following a Maximum a Posteriori approach, estimate of 
the original data is obtained by minimization of the cost {'unction given by Equation 
2.17. To provide this, a smoothing technique is devised which is based on a 
combination of Kalman Filtering (KF) and Dynamic Programming { OP t. Dynamic 
Programming is used to determine the best sequence of edges which maximizes the 
likelihood function and a detailed algorithm is presented in Appendix A. Kalman 
Filtering is used to filter within the regions. This filtering technique is preferred 
because the Kalman Filtering is the best linear optimal estimator. If the noise in the 
data is assumed to be Gaussian, the Kalman Filter gives the minimum variance 
estimate of the original scene. In particular it evaluates the conditional mean of the 
original data given the past observations (measurements). If the assumption of 
Gaussian noise is removed, the KF yields the best linear estimate. 

In the next subsection, the segmentation algorithm for a binary sequence (i.e.. 
one which assumes only two levels) is presented and the result will be extended to 
general multilevel cases. Related computer programs are included in Appendix B. 

A. BINARY SEGMENTATION USING DYNAMIC PROGRAMMING 

ALGORITHM 

Suppose that it is desired to segment a binary sequence having logic levels "1" 
and "<)" which are related to two intensity levels (gray levels) corrupted by an additive 
noise. The noise is considered as a zero mean Gaussian noise with a known standard 
deviation cr. In the general image processing extension, we can consider a logic "0" as 
a background and a logic "1" as an object. An example of the original sequence and 
the noisy sequence is given in the next chapter. 

Let Xj be the logic values of the original data of intensity levels, which assigns 
real numbers F(0), F(l) to the regions denoted by logical zeroes and ones respectively. 
Then the noisy data can be given as 

h l = F {Xj ) + W. (3.1) 

where W. — X (0 , a). We can define the signal to noise ratio (SNR) as 
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SNR = 



F( 1) - F(0) 



<? 






In this section we assume to know the levels r(U). F(l) and the noise variance t»-. 
Here, the assumption of F and <7 known is not restrictive, since various algorithms 
exist in the literature which enable one to estimate means and variance of mixtures of 
Gaussian populations. We have seen in Chapter 2 that a MRF can be arbitrarily 
assigned by the potential function V. In particular we can model spatial continuity by 
assigning hieh probabilities to smooth signals. This can be achieved by assigning the 
likelihood function in the following form: 



N, 



f (V x i x N 1 )=Pl ? (x i . 

i = 0 



, "> 

rV --H>rF(V 

i = 0 



/ ■> V 



where 



g(x, , x 2 ) = 



1 if Xj = x, 

-1 otherwise 



( 3 . 4 ) 



The parameter P models the smoothness of the original data. In recursive form 
the likelihood function can be written as: 



f k+ l( x 0’ x W"’ x k+ l) = E k( x 0- x l’-’ x k) + Ps(x k+ i- x k )-4T2~l>k+ r F ( x k+ i >l 2 (3-5) 



with C t = 0 . The estimate x is obtained by maximizing f over all possible realizations 
x.. Applying Dynamic Programming! DP) algorithm in Appendix A, the maximum of 
f k is obtained. Results of segmentation of noisy binary image is given in Chapter 4. 
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B. SMOOTHING OF PIECEWISE CONSTANT SIGNALS IN ONE 
DIMENSION USING KALMAN FILTERING 

The algorithm presented in the previous section foi binary data is suitable tor a 
small number of levels. Although it can be generalized to any number of levels, its 
complexity increases considerably and becomes unfeasable in multilevel situations. So. 
tor one and two dimensional observations which have more than two intensity levels, 
we need to search tor a ditferent approach. 

As mentioned in the first section of this chapter, a new algoiitlnn which uses the 
Kalman Filtering techniques and Dynamic Programming is developed. In applications 
filtering is concerned with the extraction of signals from noise. If the observation and 
original scene are modeled by linear stochastic models, a solution to the general 
filtering problem can be obtained by use of the Kalman Filter. 

Due to the assumed piecewise characteristics of the data, the realizations of V 
can be modeled by the state-space models 



x t+l = 



X t 4- v t 



(3.6) 



>t = 



w. 



(3.7) 



where w is Gaussian zero mean i.i.d. with standard deviation <7. and v is nonzero at 
the edges between regions only. By defining the original data X in one dimension as 

X = (X t . t e L], L = (1.2 N}. the one dimensional Markov model corresponding to 

MRF has joint density 



lnP x (.\)=p^2( x t> x t-i ) ' lnZ 

teL 

and the likelihood function to be maximized like Equation 2.16 is 



(3.S) 



f(x)= - + Pls( x t-Vt ) 



(3.9) 



teL teL 

where the first term on the right hand side of this equation comes from the noise and 
the second term is the potential defined before. 
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The Kalman Filter for the state-space model in Equation 3.6 and 3.7 is defined 
by the recursion Equation 3.11 as 



t-M 



A 

= X. 



+• K 



t+ 1 



(y t -V 



( 3 . 10 ) 



where K is the Kalman gain. When an edge is detected the gain needs to be 
reinitialized at that point. To do this, define a binary sequence 



m t = v t (3.11) 

which has a logic "1" at the edges between the regions and a logic "0" in the regions. 
Just to give an example, let's say Tj = [tj. U • 1] and T, = [t., , t 3 - 1] where T ( . T-, are 
adjacent regions and tj < t-, < t^. We can define = v ( , = 1 as shown in 

Figure 3.1 In general 



v 

t 



1 if t - 1 e T k and t s T k for some k. 
0 otherwise 



(3.12) 




Figure 3. 1 An Example for Adjacent Regions. 
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Then the estimated data can be determined from the noisy sequence y mine Equatum 
3. 10 and 3. 1 1 as 



/S 

X. 



/N 

= X. 



K. 



i ( -YV 



( 3 . 13 ) 



where 




( 3 . 1 - 4 ) 



and. 



c 



t 



1 ifm t = 1 
T j + 1 otherwise 



(3.15) 



Filtering liquation 3.13 yields the estimate of x as 



x. = 



is- 



j=l 



( 3 . 10 ) 



where t stands for the distance of t from the nearest left edge. liquation 3.10 
expresses the Kalman Filter as an averaging of the measurements. 

In the case that the edges are not known, observe the following: 

1. Given the noisy data v for any binary sequence m, the estimate x is well defined 
by Equation 3.13 - 3.15, call this estimate x(mt: 

2. The g(. , .) term in the Markov model Equation 3.S and the likelihood function 
Equation 3.9 can be interpreted as a penalty to the edges. From the 
dependence of x on the sequence m stated above it can be written. 




1 if m t = 0 
-1 ifm t = 1 



Therefore we can define the function 



( 3 . 17 ) 
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Y( n \) = 



(3. IS) 



1 il' m t = 0 
-1 ifm t = 1 

3. For any j e L let y>, xh xE m J be the sequences and y. x, x. m up to index j. for 

instance y J = {y () , Vj y. } . With this definition the likelihood function 

Equation 3.9 associated with the estimate of x satisfies the recursion 

C ) + 1 (m J + 1 )=C ) (m^) -^-{y. + 1 -x j + 1 (m ) + 1 )) 2 + Py(m j + l ) (3.19) 

In the case of the index set is L = (0,1 N} the likelihood function Equation 

3.9 for the estimate x(m) is 

f(x(m)) = (^(xfm)) (3.20) 

By these considerations the smoothing problem of the data is reduced to 
maximization of the likelihood function Equation 3.9 with respect to all binary 
sequences m. This can be done using the recursive formula Equation 3.19 and 
Dynamic Programming techniques given the details in Appendix A. 

C. SMOOTHING OF PIECEWISE CONSTANT SIGNALS IN TWO 
DIMENSIONS USING KALMAN FILTERING 

There is a problem in extending Kalman Filtering techniques from one dimension 
into two dimensions due to the lack of a causal state-space model for higher 
dimensions. Anyway, we can determine a two dimensional recursive formulation like 
Equation 3.13 - 3.15 by assuming the estimate x as the average value of the noisy data 
within the regions. 

Proceeding just like the case in one dimension, a smoothing algorithm can be 
developed to maximize the likelihood function given below: 

N N X 

£(k,t)= -^Il^k f t)-x(k f t)| 2 + py g (x ktt ,x ktt . 1 )+p£g(x ktt .x k . 1 ^ (3.21) 

k,t k,t k,t 

To do this, first assume the edges to be known, and define M = (m k } where m k 6 
(e 0 .ej.e^.e,} for each (k.t) e L. e. ( being the four combinations of edges corresponding 
to the four possibilities in clique q k t l as 

•>■> 



(3.22) 



and 

8< x u • x k-u> - * 1 

These combinations of edges is shown in Figure 3.2. 



(3.23) 




Figure 3.2 Illustration of Edges. 

Some definitions are given below referring to Figure 3.3 for each point (k,t) c L to 
obtain the smoothing algorithm: 

Definition /: r k ( is the number of points between (k.t) and the closest edge 
above. 

Definition 2: z k t is the average value of the noisy data y in the (t, . x 1) region. 
Definition 3: k k t is the number of points in the homogeneous region surrounded 
by row k. column t and the line of edges. 

Definition 4: x k t is the average of noisy data y in the region in which a., t is. 
Using the definitions given above, the filtering algorithm can be developed as 
following 



z k,t z k-l,t ~ ^k.t^k.f z k-l.f 



(3.24) 
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Figure 3.3 Illustration of Parameters of Two Dimensional Kalman Filter. 



where 



K 



k.t 




(3.25) 



and 



where 



(3. 2d) 



H 



k.t 




(3.27) 



also 
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1 



.-> } 



t 



k.t 



T k-i.t + 1 lfin k.t e 

1 otherwise 



1 '*0 



e 



V 



K 



k.t 



^'k.t-1 + T k.t ^ nl k,t e ‘ C fj ’ c 2' 

t otherwise 



(3.29) 



This algorithm is based on the fact that for any disjoint pair of regions A . B c L, the 
average of a noisy signal denoted by x on the regions A . B , A U B are related as 

x(A U B) = ajX(A) + a,x(B) ( 3.30 j 



with 3j = |A| !A U B|, a, = |B| |A U B| where |.| denotes the number of elements in 
the set. The recursion Equation 3.24 computes the average on A and Equation 3.26 
the average on A U B. As in the one dimensional case Equations 3.21 - 3.24 represent 
a well-defined mapping which associates an estimate \(M) to any field of edges M e 

v x y 

je^epe^.e,}- l • 2 . Analogous considerations as for the one dimensional case lead 
to the likelihood function as 







x N 2 (xlm(v u )]} 



k 3.3 1 ) 



where 



r 



m < v k,t } = { 



* ifv k,i = e o 

0ifv k,t = e l ’ e 2 
h * l ifv k,i = e 3 



So, the estimate x can be determined by searching over all possible sets of edges M in 
order to find x which maximizes the likelihood function Equation 3.21. Hcwe\er the 
Dynamic Programming algorithm is used to maximize the likelihood function in two 
dimensions. 
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IV. RESULTS OF SIMULATIONS 



A. ONE DIMENSIONAL SIGNALS 

We started the simulations by creating a one dimensional two level binary 
sequence which has 128 points as shown in Figure 4.1. Gaussian zero mean noise 
with different signal to noise ratios (SNRs) calculated with the formula given by 
Equation 3.2 was generated and added to the original signal by using IMSL library 
subroutines. The noisy sequences are given in Figure 4.2. 4.3, 4.4 and 4.5. Here, the 
standard deviations for the noise are 25. 50, 75 and 100 respectively. 




Figure 4.1 Original Binary Sequence. 
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Figure 4.2 Noisy Sequence with c= 15 
SNR = 4. 




Figure 4.3 Noisy Sequence with <7=50 
SNR = 2. 
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Figure 4.4 Xois\ Sequence with <7=7 5 
SNR = 1.25. 
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Figure 4.5 Noisy Sequence with <7 = 100 
SNR = 1. 

The segmentation algorithm explained in the first section of Chapter 3 which uses 
Dynamic Programming was applied to iiiter out the noisy signal with the signal to 
noise ratios given above. To improve the resuit, for each simulation the algorithm 
was run twice, from left to right and from right to left in order to provide smooth 
results. 

The outcome of simulations depend on the model parameter p that models the 
smoothness of the original data and this parameter was switched several times until we 
had the best segmentation for each case. The values of |3 were set by trial and error 
for this study. The best values for P for the given SNRs were obtained as >),34. 0.80, 
0.~4 and 0.72 respectively. Filtered signals are presented in Figures 4.0. 4.7. 4.S and 
4.9. 
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Figure 4.6 Segmented Sequence for < 7=25 and JJ = 0.S4. 
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Figure 4.7 Segmented Sequence for <7 = 50 and P = 0.50. 



s 

J 



I 

I 




Figure 4.S Segmented Sequence for < 7 = 75 and |) = 0.~4. 




Figure 4.9 Segmented Sequence for <"= 100 and P = 0.72. 

For the case of data with unknown levels we developed a different segmentation 
algorithm. The difficulty is to estimate the data levels, as well as to determine the 
segmentation. This segmentation algorithm is based on Kalman Filtering and 
Dynamic Programming. First, this method has been applied to tne one dimensional 
noisy signals shown in Figures 4.2. 4.3. 4.4 and 4.3. The noisy data was filtered 
twice (two passes). The two passes results are shown in Figures 4.10, 4.11 for <7=25 
and 30. 
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Figure 4.10 Segmented Sequence by KF with cr 
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Figure 4.11 Segmented Sequence by KF with <T = 50. 

B. TWO DIMENSIONAL IMAGES 

We applied the two dimensional filtering algorithm to the test images which 
contain three different regions given in Figure 4.12. In this test image there are four 
different intensity levels (one for the background and three for the objects). Random 
Gaussian noise has been added with standard deviation <J= 10 and (7=20 as shown in 
Figures 4.13, 4.14. 

The result of filtering is shown in Figures 4.15 and 4.16. Notice that in these 
cases the noise has been removed while preserving the edges between the regions. The 
values of the parameter P giving the best results are P=2.0 for cr= 10 and P= l.S for 
<7 — 20 showing a trend of this algorithm. The value of P depends on the SNR of the 
noisy picture. In particular p should be decreased for higher levels of noise. 
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As an application to a test image in underwater environment, we applied the 
algorithm to a 512 x 512 picture of a fish corrupted by additive noise. The noisy 
image is shown in Figure 4.17 and the filtered one is shown in figure 4. If The 
significance of it is the fact that after filtering the fish and the background present well 
distinct intensity levels. Although applications to this class of problems are still under 
investigation, the fact that the object (.the fish in this case) presents characteristics well 
distinct from the background can be used for automatic detection or recognition m an 
artificial intelligence framework. 




Figure 4.12 Original Test Image. 
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Figure 4.13 Test Image with <?= 1U. 
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Figure 4.14 Test Image with <7 = 20. 
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Figure 4.15 Segmented Image with <r= 10 and P = 2.0. 
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Figure 4. 17 linage of Fish with <T = 25. 
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Figure 4.18 Segmented Fish with a =25 and P = 0.25. 

Analogously we tested the algorithm on a 16 level checkboard shown m Figure 
4.19 for different levels of noise. In particular the values of the 16 levels are given by 
10". 50, 180, 70, 200, 120, 60, 140, 75, 175, 90, 65, 130, 55. 190, 110 from left to right 
for each line. The effects of additive noise are given in Figures 4.20 and 4.21 for <7 = 10 
and <7=20 respectively. The noise is always assumed to be Gaussian, zero mean md 
white. The application of the filtering algorithm is shown in Figures 4.22 and 4.23 
using P=l.2 and P = 0.7 respectively. As expected the algorithm filters within the 
regions while preserving the sharp separation between adjacent regions. Also shown in 
Figure 4.22 are the edges between regions, as detected by the algorithm which shows 
the reinitializations of the Kalman Filter gains. 
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Figure 4.19 Test Image with 16 DitTerent Regions. 
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Figure 4.20 Test Image with 10. 
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Figure 4.21 Test Image with <T= 20. 
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Figure 4.22 



Segmented Image with c= 10 and 1.2. 
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Figure 4.23 Segmented Image with <7 = 20 and p = 0.7. 
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V. CONCLUSIONS 



An algorithm to segment noisy data in one dimension and two dimensions has 
been presented, where the segments are piecewise constant. The data are described by 
state-space models. The particular feature of the algorithm is to search lor the best 
sequence of edges in order to maximize the likelihood function. The algorithm has 
emphasized piecewise constant data. bur. it can be extended for the data with legions 
characterized by textures to which we associate different Autoregressive (AR) models. 
Typical applications are not only image segmentation but also segmentation in speech 
analysis, such as phenomenon separation. 
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APPENDIX A 

DYNAMIC PROGRAMMING ALGORITHM 



The Dynamic Programming is an algorithm used when the solution to a problem 
may be considered as the result of a sequence of decisions. An optimal sequence of 
decisions will maximize the given function, in the case that is used in this study, the 
likelihood function f k . 

In Dynamic Programming an optimal sequence of decisions is obtained by 
making explicit appeal to the Principle of Optimality. In the cases when this principle 
can be applied, an optimal sequence of decisions has the property that whatever the 
initial state and decision are. the rest of the decisions must make up an optimal 
decision sequence according to the state resulting from the first decision. Standard 
Dynamic Programming techniques are introduced by Bellman. [Ref. 9| 

The mathematical equations and the algorithm is given below for the case that 
was mentioned in the second section of Chapter 3. By defining the likelihood function 
as 



f k (x 0 .Xi,...,x k )=pVg(x r x. i> - — Ily.-F(x )| 2 
j = 0 2<7 2 j = 0 



or in recursive form 



f k+i( x o- x i x k+i)=Mv x i x k) +Ae k+i( x k+i’ x k) 



where 



(AT) 



(A. 2) 



A W x k+i’ x k )=Pg(\+i>\)- 



“~l>k+r F ( x k + i)l 2 
2<t 2 



1A.3) 



and setting the initial condition = 0 and also assuming that x k can be only logic 

"0" or "1", we can determine the best sequence {x k , k = 0,...,N} which maximizes 
(v s (\ u ,...,Xy). For this purpose define 
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< -V4) 



V°) = max f k (\ 0 x L .,.U) 



Vf 



^ X 0 X k-1 1 



and 

JA 1 ) = max f k (x Q x k _j. 1) (A. 5) 

( v 'o W 

We can represent the problem in the form of a graph as in figure A.l. with N stages 
where the nodes at the k-th stage represent the state x k (either "0" or "l"), and the 
branches represent the gain associated to the transition from one state to the next. In 
this way J k (0) and J k ( 1 ) represent the maximum of all possible gains up to x k = 0 and 
x k = 1 respectively. These quantities can be recursively updated as 

J k+ j(0) = maxj J k (0)+ AC k + ^ (0,0 ) , J k ( 1 )-r AC k + ^ 1,0)] t A. 6) 

and 

J k+1 (l) = max!J k (0) + AC k + 1 (0.1) ,.I k (l) + Ae k+1 (l,l)} < A.7> 

with AC k as in Equation A. 3. Also we can keep track of the branches which yield the 
maximum likelihood by defining Pointer k (0), and Pointer^ 1) at each stage k. The 

sequence x maximizing C\;(x 0 x^-) is therefore obtained by backtracking as in the 

following: 

Let x^- be such that Jn^x^) = max[J>q(0) . 1)] 

for k = N - 1 to 0 

\ = Pointer k + 1 (x k+ ,) 

end for 

end 
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k=0 

logic * 0 



k=1 27 



START 




Figure A. I Choice of Best Sequence of Decisions. 
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APPENDIX B 

COMPUTER PROGRAMS 






* 

* 

k 

k 

k 

k 

k 



PROGRAM BINARY . DAT 

PURPOSE To generate 128 point binary sequence. 
REQUIRED IMSL ROUTINES None. 

IMPLEMENTED BY Lt. Kani HACIPASAOGLU Jan. 1987 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkx:kkkkkkkkkkkkkkk7<;kk7r7S7tkk 



INTEGER IN ITT , FI MITT , F( 128) , TERM,SETEUF,K 
REAL DWINDO , MOVEA , DRAWA 

OPEN(UNIT=l,NAME=' BINARY. DAT 1 ,ACCESS= 1 DIRECT 1 , STATUS= ' NEW ' , 
* RECL=32,MAXREC=123) 

C 

DO 10 K=1 , 128 

IF (K .GE. 0 .AMD. K .LT. 30) THEM 
F(K)=100 

ELSE IF(K .GT. 70 .AMD. K .LT. 90) THEN 
F(K)=100 

ELSE' IF (K .GT. 120 .AMD. K .LE. 128) THEN 
F(K)=100 
ELSE 

F(K)=200 
END IF 

WRITE (1 1 K) (F (K) ) 

10 CONTINUE 

CLOSE (UMIT=1) 

C 

CALL INITT(480) 

CALL TERM (2,1) 

CALL SETBUF ( 2 ) 

CALL DWIMDO( -50. 0,300. 0,-50. 0,400.0) 

CALL MOVEA (-50.0, -50.0) 

CALL DRAWA(300 .0 , -50 .0 ) 

CALL DRAWA (300.0,400.0) 

CALL DRAWA (- 50. 0,400.0) 

CALL DRAWA (-50.0,-50.0) 

CALL MOVEA (0 .0,1000.0) 

CALL DRAWA (0 . 0 , 0 . 0 ) 

CALL DRAWA(128. 0,0.0) 

CALL MOVEA ( 0 . 0 , 0 . 0 ) 

DO 35 K=1 , 128 
X=K 
Y=F(K) 

CALL DRAWA (X,Y) 

35 CONTINUE 

CALL FINITT(0 ,0 ) 

STOP 

END 
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* k 

* PROGRAM B I NARY 2 . DAT * 

k k 

* PURPOSE To generate one dimensional 128 point binary * 

sequence v/ith zero mean Gaussian" noise use 
by Dynamic Programming algorithm to segmen 

* a noisy signal. * 

k k 

* REQUIRED IMSL ROUTINES GGNHL * 

k k 

* IMPLEMENTED BY Lt. Kani HACIFASACGLU Feb. 1937 * 

k -k 

k-k'kkkkkk'kkk'kkkkkk'kkkkkkkkkkkk.-kkkk-kkkkk'kkkkkkkkkkkkkkkk kkkkk kk 



INTEGER INIIT,FINITT,F(12S) , TERII , SETBUF , I , MR , L 
REAL DWIHDO , MOVEA , DRAWA , GGNHL , R , S ( 128 ) 

DOUBLE FRECISION DSEED 
DSEED=65471 

OPEN (UNIT=2,NAME=' BINARY2.DAT 1 , ACCESS= 1 DIRECT ' , STATUS* 1 NEW 1 , 
* RECL=32,MAXREC=128) 

C 

DO 10 1=1,128 

I F ( I .GE. 0 .AMD. I .LT. 30) THEN 
F ( I )=1O0 

ELSE IF ( I .GT. 70 .AND. I .LT. 90) THEN 
F ( I ) =100 

ELSE IF ( I .GT. 120 .AND. I .LE.128) THEN 
F ( I ) =100 
ELSE 

F ( I )=200 
END IF 
NR=1 

CALL GGNHL (DSEED, NR, R) 

SIGMA=7 5 . 0 
R=SIGHA*R 

C R : NOISE FUNCTION WHICH WAS GENERATED 
S(I)=F(I)+R 
WRITE (2 1 I) (S(I)) 

10 CONTINUE 

CLOSE(UNIT=2) 

C 

CALL INITT(480) 

CALL TERM (2,1) 

TAT T qFTRTIF ( ?\ 

CALL DWINDO (-50.0,300.0,-50.0, 400 . 0 ) 

CALL MOVEA(- 50.0, -50.0) 

CALL DRAWA ( 300.0, -50.0) 

CALL DRAWA ( 300. 0,400.0) 

CALL DRAWA (-50. 0,400.0) 

CALL DRAWA (-50.0, -50.0) 

CALL H0VEA(0. 0,1000.0) 

CALL DRAWA ( 0 .0,0.0) 

CALL DRAWA(l28. 0,0.0) 

CALL M0VEA(0. 0,0.0) 

DO 35 1=1,128 
X=I 

Y=S(I) 

CALL DRAWA (X,Y) 

35 CONTINUE 

CALL FINITT(0 , 0 ) 

STOP 

END 
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* PROGRAM 0MUR2 - 

k k 

PURPOSE To segment a noisy binary sequence using 
Dynamic Programming algorithm. 

* REQUIRED IMSL ROUTINES None. * 

* IMPLEMENTED BY Lt. Kani HACIPASAOGLU March 19S7 

* k 
^kkkkkkkkkkkkkkkkkkkkkkkkkkkrkkkkkkkkkkkkkkkkkkkkkkk7tkkkkkkkkkk-kk 



INTEGER H , K, FO , FI ,OUT , PTR( 128,0:1) ,F(123) 

INTEGER INITT , FINITT , TERM , SETBUF , J , M 

REAL BETA , S IGMA , LNEWO , LNEW1 , MOVEA , L0 , LI , OUTS (1:123) , SUM ,T" 
REAL DELV (1:123,0:1,0:1), DWIMDO , DRAW A , S ( 128 ) , X , Y , S3E , SEE , T 
S IGMA=1 50 . 0 
READ ( * , 1 ) BETA 
1 FORMAT (F4.1) 

F0=10Q 

F1=2Q0 

OPEN(UNIT=2 ,MAME= 1 BINARY2 . DAT 1 , ACCESS= 1 DIRECT 1 , STATUS= 1 OLD 1 , 

* RECL=32 , MAXP.EC=1 28 ) 

DO 33 M=1 .123 

READ (2 'll) ( S (N) ) 

WRITE ( 6 , *) S(N),N 
33 CONTINUE 

CLOSE (UNIT=2) 

K=1 

L0=BETA-( ((S(K)-100)**2)/(2*(SIGMA**2) )) 

Ll=-EETA-(( (S(K)-200)* ;1: 2)/(2 :r (SIGMA* *2 ) ) ) 

DO 15 H=1 , 123 
K=H- 1 

DELV(K+l,0,0)=BETA-( ((S(K+1)-F0 )**2)/(2*( SIC-MA**2))) 
DELV(K+1, 1,0 =-BETA-(( ( S (K+l ) -F0)**2)/ ( 2*(SIGMA**2 ) ) ) 
DELV(K+l,0,l)=-BETA-( ( (S (K+l ) -FI ) **2) / ( 2* (SIGNA**2 ) >) 
DELV(K+l,l,l)=BETA-( ( (S (K+l) -FI )**2 )/ (2* (SIGHA**2) )) 

IF ( (L0+DELV(K+1 ,0,0)) .LT. (L1+DELV(K+1, 1,0) ) ) THEM 
LNEW0=L1+DELV(K+1 ,1,0) 

PT R ( M , 0 ) = 1 
ELSE 

LNEW0=L0 +DELV(K+1 ,0,0) 

PTR(M,0)=0 
END IF 

IF( (L0+DELV(K+1 ,0,1)) .LT. (L1+DELV(K+1 , 1 , 1) ) ) THEN 
LNEW1=L1+DELV(K+1 ,1,1) 

PTR(M,1)=1 

ELSE 

LNEW1=L0+DELV (K+l ,0,1) 

PTR(M, 1)=0 
END IF 

L0=LNEW0 
L1=LNEW1 
15 CONTINUE 
C TRACE BACK PROCEDURE 

OPEN(UNIT=3 ,NAME= 1 OUTS . DAT 1 ,ACCESS= 1 DIRECT 1 , STATUS= 1 NEW 1 , 

* RECL=32 ,MAXREC=128) 

IF (L0 .LT. LI) THEN 
J=1 
ELSE 
J=0 
END IF 
OUT=J 

DO 44 K=127 ,0,-1 
IF (OUT .EO. 0) THEN 
OUTS (K+l) =100 
ELSE 

OUTS(K+1)=200 
END IF 
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0UT=PTR(K+1 , OUT ) 

WRITE (6 , *) OUTS (K+l ) , K+l 
44 COMTIMUE 

DO 55 K=1 , 128 
WRITE ( 3 1 K ) (OUTS (K) ) 

55 CONTINUE 

CLOSE (UMIT=3) 

C 

CALL INITT (480 ) 

CALL TERM (2,1) 

CALL SETBUF ( 2 ) 

CALL DWIMDO (-50.0,300.0,-50.0,400.0) 
CALL IlOVEA(- 50.0, -50.0) 

CALL DRAWA( 300 . 0 , -50 . 0 ) 

CALL DRAWA( 300. 0,400.0) 

CALL DRAWA(-50. 0,400.0) 

CALL DRAWA(- 50.0, -50.0) 

CALL HOVEA( 0.0, 1000.0) 

CALL DRAWA(0 .0,0.0) 

CALL DRAWA( 150.0,0.0) 

CALL MOVEA(0. 0,0.0) 

DO 3 K=1 , 128 
X = K 

Y=OUTS (K) 

CALL DRAWA(X,Y) 

3 COMTIMUE 

CALL FINITT (0,0) 

STOP 

END 
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* IT 

* PROGRAM TEST . DAT * 

* TT 

PURPOSE To generate a test image which has * 

three different intensity levels. A 

* 1 ~k 

* REQUIRED IMSL ROUTINES None. * 

A X 



* IMPLEMENTED BY Lt. Kani HACIPASAOGLU May 1937 * 

* rr 

BYTE A( 123 , 123) 

INTEGER R,XC1, YC1,XC2,YC2, FI (1:128,1: 128) ,F2(1: 123,1 : 123) 

R=20 
XC1=60 
YC1=60 
YC2=60 
XC2=90 

OPEN ( UNI T=5 , NAME-” 1 TEST . DAT 1 , ACCESS= 1 DIRECT 1 , STATUS= 1 NEW 1 , 
* RECL=32,HAXREC=128) 

C 

DO 10 1=1,123 
DO 20 J=1 , 1 23 

F1(I,J)=(((I-XC1)**2+(J-YC1)**2)-(R**2) ) 

F2(I, J) = ( ( (I-XC2)* J: 2+(J-YC2)**2)-(R :r *2) ) 

IF( ( FI (I , J) .LE. 0) .AND. (F2(I,J) .GT. 0)) THEN 
A(I, J) = 50 

ELSE IF((F1(I,J) .GT. 0) .AND. (F2(I,J) .LE. 0)) THEN 
A(I , J) = 150 

ELSE IF( (FI ( I , J) .LE. 0) .AND. (F2(I,J) .LE. 0)) THEN 
A ( I , J ) =100 
ELSE 

A(I, J)=200 
END IF 
20 CONTINUE 

WRITE ( 5 1 1 ) ( A( I , J) , J=1 , 128 ) 

10 CONTINUE 

CLOSE(UNIT=5) 

STOP 

END 
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PROGRAM TESTIM.DAT 

PURFOSE To generate a test image with zero mean 
Gaussian noise to be segmented by 
Kalman Filter. 

REQUIRED IMSL ROUTINES GGNHL 

IMPLEMENTED BY Lt. Kani HACIPASAOGLU May 1937 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkx 



BYTE A(123, 123) 

CHARACTERS TESTIM 
INTEGER XC1 , YC1 , XC2 , YC2 , FI ( 128 , 128) 

INTEGER NR, Rr , SIGMA, CR,F2 (123, 128 ) ,AA( 123 ,128) 

REAL GGNML , R 
DOUBLE PRECISION DSEED 
DSEED=6547 1 
WRITE (*,100) 

100 FORMAT ( 'ENTER IMAGE DATA NAME 1 ) 

READ(*, 101) 

101 FORMAT (A50) 

C 

CR=20 

XC1=60 

YC1=60 

YC2=60 

XC2=90 

C 

OPEM(UNIT=l ,NAME=' TESTIM. DAT 1 ,ACCESS= ' DIRECT ' ,STATUS= 1 NEW 1 , 
* RECL=32,MAXREC=128) 

C 

WRITE (*, 222 ) 

222 FORMAT ( 'ENTER SIGMA' ) 

READ (*,333) SIGMA 
333 FORMAT ( 12 ) 

DO 10 1=1,128 
DO 20 J=1 , 128 

FI (I , J)=( ( (I-XC1)**2+(J-YC1)**2)-(CR**2) ) 

F2(I, J)=( ((I-XC2)**2+(J-YC2)**2)-(CR**2) ) 

IF ( ( FI ( I , J ) .LE. 0) .AND. (F2(I,J) .GT. 0)) THEN 
A( I , J) = 50 

ELSE IF( (FI ( I , J) .GT. 0) .AND. (F2(I,J) .LE. 0)) THEN 
A( I , J)=l 50 

ELSE IF ( (FI ( I , J) .LE. 0) .AND. (F2(I,J) .LE. 0)) THEM 
A( I , J)=100 
ELSE 

A( I , J)=200 
END IF 
NR=1 

CALL GGNML (DSEED.NR,R) 

Rr=SIGMA* jnint (R) 

AA( I , J )=A( I , J)+Rr 
I F ( AA ( I , J ) .LT. -128) THEN 
AA(I , J)=AA(I , J)+256 
ELSE I F ( AA ( I , J ) .GT. 127) THEN 
AA(I , J)=AA( I , J) -256 
ELSE 

AA( I , J)=AA( I , J ) 

END IF 

A( I , J )=AA(I , J ) 

20 CONTINUE 

WRITE ( 1 1 1 ) (A(I , J) , J=1 , 128) 

10 CONTINUE 

CLOSE (UMIT=1) 

STOP 

END 
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k -4- 

* PROGRAM FINAL . DAT 

k 7 ? 

PURPOSE To generate a test imaye which has 16 ** 

* different regions. 

k 

* REQUIRED IMSL ROUTINES None. 

x k 

* IMPLEMENTED Ex Lt. Kani HACIPA3AOGLU June 1987 * 

k r. 

kkkkkkkkkkkkkkkkkkkkkkkkkkkk kkkkkkkkkkkkkkkkkk kkkkkkxkkkkk k7tx 

BYTE A( 123 , 128) 

C 

0FEI!(UMIT=1 , MAME= 1 - INAL . DAT ' , ACCE3S= ' DIRECT 1 , STATUS= ' NEW 1 , 

* RECL=32 ,I1AMREC=128 ) 

DO 10 K=1 . 32 
DO 20 L=1 , 32 

A(K,L)=100 
A(K,L+32)=50 
A( K , L+64 )=130 
A(K,L+96;=70 
C 

A(K+32 , L)=200 
A(K+32,L+32)=120 
A(Kf32,L+64)=60 
A(K+32,L+96)=140 
C 

A(K+64,L)=75 
A(K+64,L+32)=175 
A< K+64 , L+64) =90 
A(K+64 , Li-96 )=65 
C 

A(K+96 ,L)=130 
A(K+96 ,L+32)=55 
A(K+96,L+64)=190 
A(K+96,L+96)=110 
20 CONTINUE 

10 CONTINUE 

DO 30 1=1,128 

WRITE ( 1 ' I ) (A( I , J) , J=1 , 123) 

30 CONTINUE 

CLOSE (UMIT=1 ) 

STOP 

END 
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* * 

* PROGRAM NFINAL.DAT 

* x 
A 
;*: 

A 

A 

* 

A 
A 

A A 

a a aaa a a aa*^* a a a a kkkkkkkkkkk a A^r k kkkkkkkk kkk kkkk k a* a** a a* a a a a a 



PURPOSE To generate a test image which has 16 

different regions with zero mean Gaussian 
noise to be segmented by Kalman Filter. 

REQUIRED IMSL ROUTINES GGMHL 

IMPLEMENTED BY Lt. Kani HACIFASAOGLU June 1937 



BYTE A(123 , 128) 

CHARACTER*5Q MFIMAL 

INTEGER MR, Rr , SIGMA, AA( 128, 128) 

REAL GGMHL , R 
DOUBLE PRECISION DSEED 
DSEED=6547 1 
WRITE ( * , 100 ) 

100 FORMAT < 1 ENTER IMAGE DATA NAME') 

READ(*,101) 

101 FORMAT (A50) 

OPEN ( UNIT=1 ,MAME=' NFINAL3.DAT' ,ACCESS= 1 DIRECT 1 , STATUS= 1 MEW 1 , 
* RECL=32 ,MAXREC=128) 

WRITE ( * , 222 ) 

222 FORMAT ( ' ENTER SIGMA 1 ) 

READ(*,333) SIGMA 
333 FORMAT (12) 

DO 10 K=1 , 32 
DO 20 L=1 , 32 

A(K,L)=100 
A(K , L+32)=50 
A(K,L+64)=180 
A(K, L+96 )=70 

A(K+32,L)=200 
A(K+32,L+32)=120 
A(K+32,L+64)=60 
A(K+32 , L+96 )=140 

A(K+64,L)=75 
A(K+64,L+32)=175 
A(K+64,L+64 =90 
A(K+64,L+96)=65 

A(K+96 ,L)=130 
A(K+96,L+32)=55 
A(K+96,L+64)=190 
A(K+96 , L+96) =110 
20 CONTINUE 

10 CONTINUE 
NR=1 

DO 80 1=1,128 

DO 90 J=1 , 128 

CALL GGNML ( DSEED , NR , R ) 

Rr=SIGMA*jnint (R) 

AA( I , J )=A ( I , J )+Rr 
I F ( AA ( I , J ) .LT. -128) THEN 
AA(I , J)=AA(I , J)+256 
ELSE I F ( AA ( I , J ) .GT. 127) THEN 
AA( I , J )=AA( I , J) -256 
ELSE 

AA ( I , J ) =AA ( I , J ) 

END IF 

A ( I , J ) =AA ( I , J ) 

90 CONTINUE 

WRITE ( 1 ' I ) (A( I , J) , J=1 , 128) 
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CONTINUE 
CLOSE (UNIT=1 ) 
STOP 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



777 

732 

773 

c 

779 

730 

790 

791 



792 
310 

793 
811 

794 
c 



6 

5 



c 

c 

555 

666 

c 

c 



^ x 

* PROGRAM CRISTI2D * 



* 
; k 
k 
k 
* 



PURPOSE To segment a two dimensional noisy image 
using'Markov Random Field, Dynamic 
Programming and Kalman Filtering techniques, 

REQUIRED IMSL ROUTINES None. 

IMPLEMENTED EY Lt. Kani HACIPASAOGLU Aug. 1987 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk^k 



common ny, y( 512 , 512 ), tau( 512 ), xtop( 512) , iedge , ymax 
integer start, dir, cycle 

character*2Q image, imout 
byte s (512 , 512) , v(512,512) 
write!*, 111 ) 

format (' Enter Input Image File Names') 
read !*, 778) image 
write!*, 782) 

format(' Enter Output Image File Name: 1 ) 
read!*, 778) imout 
format(a20) 

open (2, name=image, access= 1 direct 1 , status =l old' ) 
get data from file 
write (*, 779) 

format(' Enter Image Size 1 ) 
read(*, 780) npoints 
format ( i3 ) 
write(*, 555) 
read (* 666) sigma, beta 

write (*, 790) 

format( 1 see edges? (yes=l)') 
read (*, 791) iedge 
format( il ) 
write(*, 810) 
write ( * , 792) 

format(' Enter scale factor (1.0=no scale)') 
format!' output yout=yth+scale* (y-yth) ' ) 
read(*, 793) scale 
format! flO .2) 
write (* , 811 ) 

format!' 1: yth=(yav+ymax)/2 ; 2 :yth=0 . 9ymax ( - 3: enter yth; ??') 
read(*, 791) iyth 
if (iyth.eq.3) then 
write!*, 794) 
format!' Enter yth:') 
read (*, 793) yth 
endif 

do 5 k=l, npoints 

read(2’k) (S(k,j), j =1 , npoints) 
do 6 j=l, npoints 
temp=s(k,i) 
if ( temp . It . 0 . 0 ) then 
temp=temp+256 
endif 

y(k, j)=temp 
continue 
continue 
close (unit=2) 



format!' ENTER: sigma, beta') 
format ( 2f 10 .4) 

*** main loop *** 
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v/rite(*, 667) 

667 format(‘ Start Computing . ..') 

ymax=0 . 0 
yav=0 .0 
ncount=0 
nx=npomts 
ny=npoints 

gam=l .0/ (2 . 0*sigma**2) 
c 

c *** initialize 

c 

do 30 j=l ,ny 
tau( j )=Q . 0 
xtoo( j )=0 . 0 
30 continue 

c 

nx=nx- 1 

c *** main loop *** 

do 40 i=l,nx 
s tar t=l 
dir=+l 
cycle=l 

call line(i / start, dir, cycle, gam, beta, rowav) 

c 

start=ny 

dir=-l 

cycle=2 

call line(i,start, dir, cycle, gam, beta, rowav) 

40 continue 

c 

c *** compute average of output image 

yav=(ncount*yav+ rowav)/ (ncount+1; 
ncount=ncount+l 
c 

write(*, 668) 

668 format(* ... End Computing. Now transfering data to disk 
c 

c *** output data *** 

c 

if (iyth.eq.2) then 
yth=ymax*0 . 9 
endif 

if (iyth.eq.l) then 
yth=(yav+ymax)/2 . 0 
endif 
c 

write ( * , 800) 

300 format (' yave , ymax, yth =') 

write(*, 801) yav, ymax, yth 
801 format(3(x, fl0.4)) 



do 60 k=l, npoints 
do 61 j=l, npoints 

ys=yth+scale*(y(k, j ) -yth) 
c 

if (ys.qt. 255.0) then 
ys=255.0 
enaif 

if (ys.le.0.0) then 
ys=0 . 0 
enaif 
c 

itemp=jnint (ys) 
if (itemp. gt . 127) then 
itemp=itemp-256 
endif 

v(k, j )=itemp 
61 continue 

60 continue 
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c 



c 



50 

c 



c 



c 

c 

c 



90 

c 

c 



999 

c 

c 

c 

c 



c 



c 



110 

c 



c 

c 



open (unit=3, name=imout, access= ' direct 1 , 
status='new' , recl=123, maxrec=5l2) 

do 50 i=l npoints 
write(3 , j) (v(j,k), k=l, npoints) 
continue 
close (unit=3) 



subroutine line(i, start, dir, cycle, gam, beta,rowav) 
integer start, air, cycle 

common ny, y(512,512), tau(512), xtop ( 512 ) , iedge , ymax 
real xhat(4.bl2), lambda(4 , 512 ) 

real ttau(4), txtop(4), tlambda(4), txhat(4), te(4) 

real newtau(4,512) , newxtop(4, 512) 

real e (4 , 512 ) 

integer pointer (4 , 512) 

*** scan from START by DIR (=+l, -1) 
j=start 

*** initialize first column 
do 90 k=l ,4 

xhat(k, start )=y ( i, start) 
lambaa (k, start)=0 . 0 
e(k,start)=0.0 
continue 



rowav=0 .0 
ncount=0 
i=i+dir 

*** state=(left edge, upper edge) 

*** 0=no eage, l=eage 

*** state (0,0) 
do 110 k=l A 
ttau(k)=tau( j )+l 

txtcp(k)=xtop(i )+(l .0/ttau(k) )*(y(i, j )-xtop(j )) 
t lambda (k)=lambda(k. 3 -dir )+ttau(k) 
z=txtop(k) -xhat( k, j -dir ) 
txhat(k)=xhat(k, j-dir)+( ttau(k)/ t lambda (k) )*z 



if (k.eq.l) then 
emin=te (k) 
endif 

if (te(k) .le.emin) then 
im=k 

emin=te(k) 

endif 

continue 

newtau(4 , j )=txtop( im) 
lambda (4 , 3 )=tlambda ( im) 
pointer(4, j ) = im 
newxtop(4, 3 )=txtop(im) 
xhat(4 , 3 )=txhat( im) 
e (4 , 3 )=te ( im) 

*** state (0,1) *** 
do 115 k=l ,4 

tlambda(k)=lambda(k, j-dir)+l 
z=y(i, j)-xhat(k, j-dir) 
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c 



c 



115 

c 



txhat(k)=xhat(k, j-dir)+(l .0/tlambda(k) )*z 



i = < y ( i , j+dir)-txhat(k ) ) k *2 
2=|7ji,]):txhat(k);^2_ 



e 

e2 

e3=(y(i+I'j ) - txhat(k ) ) k *2 
te(k)=e(k, 3 - dir )+gam : * ; (el+e2+e3) / 3 



if (k.eq.l) then 

emin=te(k) 

endif 

if (te(k) .le.emin) then 
im=k 

emin=te (k) 
endif 

continue 



newtau( 1 , j )=1 
newxtop(l , n )=y( i, j ) 
lambda ( 1 , j )=t lambda (im) 
pointer ( 1 , 3 )=im 
xhat ( 1 , 3 )=txhat(im) 
e(l J)=te(im) 



*** state (1,0) 
do 120 k=l ,4 
ttau(k)=tau( j )+l 

txtop(k)=xtop( j ) + (1 .0/ ttau(k) )*(y (i, j ) -xtcp( j ) ) 

el=(y(i, j+dir)-txtop(k) )**2 
e2=(v( i , 3 ) -txtop(k) )^*2 
e3=(y( i+1 , j ) - txtop(k) )**2 
te(k)=e(k,]-dir) + gam*(el+e2+e3)/3 

if (k.eq.l) then 

emin=te(k) 

endif 

if ( te (k) . le . emin) then 
im=k 

emin=te(k) 

endif 

continue 



newtau(2, j )=ttau(im) 
nev;xtop(2 , j )=txtcp( im) 
lambda (2 , j ; = tau( j ) 
pointer(2 , j )=im 
xhat (2 , 3 )=txtop( im) 
e(2 ,3 ) = te(im) 



*** state (1,1) 
nev/tau(3 ,3 )=1 
newxtop(3, j)=y(i,j) 
lambda(3. j]=l 
xhat (3 , j )=y(i, j ) 
e i=fyUd + ai r )-y(i, ; i))**2 
e3=(y(i+l,3)-y(i,3)) K ^2 
z2=gam*(el+e3 



if (e(4, j -dir; . It . e ( 1 , 3 -dir ) ) then 
e(3 ,3 )=z2+e(4, j-dir)+beta 
pointer(3 , j )=4 

else 



e (3 , j )=z2+e ( 1 , 3 -dir )+beta 
pointer (3 , j )=1 

endif 



if ( j .gt. (1-dir) .and. j . It. (ny-dir) ) then 
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goto 999 
endif 
c 

c *** backtrack *** 

c 

emin=e ( 1 , j ) 

do 150 k=l , 4 

if (emin.le.e(k, j) ) then 
emin=e ( k, j ) 
im=k 
endif 

150 continue 

c 

3S3 y ( i , j )=xhat( im , j ) 

w rite(*, 111) im 
11 format(i3) 

if (cycle. eq. 2) then 

xtop(i )=newxtop(im. j ) 
tau(j }=newtau( im, j ) 

compute max intensity 
if (y(i, j ) .qt.ymax) then 
, ymax=y(i / j) 
endif 



*** row average 

rowav= ( ncoun t*rowav+y ( i , j ) ) / (ncount+1 ) 
ncount=ncount^l 

mark the edge with a black entry 

if ( im . ne . 4 . and. iedge . eq. 1 ) then 
y(i, j)=0.0 
endif 
endif 
c 

im=pointer(im, j ) 
j=j-dir 

if ( j . gt . (1-dir ) . and. j . It . (ny-dir) ) then 
goto 883 

endif 

c 

return 

end 
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